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Abstract: Despite the unparalleled diversity of venomous snakes in Australia, research has 
concentrated on a handful of medically significant species and even of these very few 
toxins have been fully sequenced. In this study, venom gland transcriptomes were 
sequenced from eleven species of small Australian elapid snakes, from eleven genera, 
spanning a broad phylogenetic range. The particularly large number of sequences obtained 
for three-finger toxin (3FTx) peptides allowed for robust reconstructions of their dynamic 
molecular evolutionary histories. We demonstrated that each species preferentially 
favoured different types of a-neurotoxic 3FTx, probably as a result of differing feeding 
ecologies. The three forms of a-neurotoxin [Type I (also known as (aka): short-chain), 
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Type II (aka: long-chain) and Type III] not only adopted differential rates of evolution, but 
have also conserved a diversity of residues, presumably to potentiate prey-specific toxicity. 
Despite these differences, the different a-neurotoxin types were shown to accumulate 
mutations in similar regions of the protein, largely in the loops and structurally 
unimportant regions, highlighting the significant role of focal mutagenesis. We theorize 
that this phenomenon not only affects toxin potency or specificity, but also generates 
necessary variation for preventing/delaying prey animals from acquiring venom-resistance. 
This study also recovered the first full-length sequences for multimeric phospholipase A 2 
(PLA 2 ) 'taipoxin/paradoxin' subunits from non-Oxyuranus species, confirming the early 
recruitment of this extremely potent neurotoxin complex to the venom arsenal of 
Australian elapid snakes. We also recovered the first natriuretic peptides from an elapid 
that lack the derived C-terminal tail and resemble the plesiotypic form (ancestral character 
state) found in viper venoms. This provides supporting evidence for a single early 
recruitment of natriuretic peptides into snake venoms. Novel forms of kunitz and waprin 
peptides were recovered, including dual domain kunitz-kunitz precursors and the first 
kunitz -waprin hybrid precursors from elapid snakes. The novel sequences recovered in this 
study reveal that the huge diversity of unstudied venomous Australian snakes are of 
considerable interest not only for the investigation of venom and whole organism evolution 
but also represent an untapped bioresource in the search for novel compounds for use in 
drug design and development. 

Keywords: venom; evolution; phylogeny; elapid; Australia; molecular evolution; 
Darwinian selection; toxin phylogenies 



1. Introduction 

Snake venoms are cocktails of toxins which have evolved from ordinary body proteins [1] to 
rapidly disrupt key physiological processes in prey animals. Since venom is energetically expensive to 
synthesize [2], an ideal venom-component would be effective even at lower concentrations. Therefore, 
target specificity of toxins is of paramount importance [3]. Since these toxins have evolved over 
millions of years of evolutionary time to rapidly and systematically breakdown prey homeostasis, they 
are invaluable as investigational ligands in elucidating physiological pathways or as lead compounds 
in drug design and therapeutics [4-8]. 

Australia is the stronghold of one of the world's most medically significant families of venomous 
snakes — the front-fanged clade Elapidae. Elapid snakes include many of the world's most infamous 
venomous snakes: the cobras of Asia and Africa; the mambas of Africa; the coral snakes of Asia and 
the Americas; the sea snakes and all of Australia's medically significant venomous land snakes. The 
Australian continent is home to at least 130 (including sea snakes) of the world's 320+ species of 
elapid snake [9]. Despite this tremendous diversity, to date, the vast majority of toxinological research 
conducted on the venoms of Australian snakes has focused on just five of Australia's 26 genera of 
terrestrial elapid snakes. Even within the five most studied genera, typically only one or two species 
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per genus has received a significant amount of attention. These five genera {Acanthophis; Notechis; 
Pseudechis; Pseudonaja; and Oxyuranus) are considered the most medically significant of Australia's 
venomous snakes [10], where "medical significance" is defined as the "danger posed to a human 
through bite". A further three genera {Austrelaps, Hoplocephalus and Tropidechis) have received a 
moderate amount of research attention, while the remaining 1 8 genera of terrestrial elapid snakes have 
historically been almost completely neglected by toxinologists. As more species have been 
investigated with novel methods, our knowledge of toxin evolution and structure-function relationships 
of these toxin types has increased. An initial investigation of the venoms of some small Australian 
elapids has shown them to be equally as complex as those of their larger, better-investigated 
cousins [11]. For this reason the small elapid fauna of Australia may be viewed as a rich and untapped 
bioresource, despite the fact that bites from many of these snakes is far from being 
"medically significant". 



Figure 1. BEAST maximum credibility ultrametric tree for in-group taxa [12]. Node 
values indicate 95% highest posterior distributions for calibration points. Posterior 
probability support values are shown for each node. Species included in this study are 
indicated in red. 
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A wide variety of toxin types have been previously sequenced from Australian elapids, including 
cysteine-rich secretory proteins (CRiSP), toxin homologues of coagulation factors Va (fVaTx) and Xa 
(fXaTx), kunitz and waprin peptides, lectin, natriuretic peptides, phospholipase A 2 (PLA 2 ), 
metalloproteinases (SVMP), and three-finger toxins (3FTx) [13,14]. However, partly due to the limited 
range of taxa sampled, our current understanding of the true distribution and diversity of 
venom-components in these snakes remains far from complete. Thus many toxin types that are now 
considered "unique" to certain species or clades may only appear to be so because of the sampling bias 
towards medically significant large Australian Elapids. 

Therefore, we investigated a wide taxonomical and ecological diversity of the under-studied 
Australian elapid snakes (Figure 1) through random sequencing of cDNA libraries; a method known to 
be efficient for biodiscovery and understanding the biochemical composition of snake venoms 
(c.f. [15-26]). The results of this study not only contribute to our understanding of the molecular 
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evolution of Australian elapid snake venom, particularly the influence of feeding ecology on venom 
composition, but will also constitute a platform for biodiscovery. 

2. Results and Discussion 

Random sequencing recovered a myriad of toxin types, previously only known from the other well 
studied species (Table 1). All venom gland transcriptomes contained sequences of multiple toxin types. 
Large globular proteins, such as acetylcholinesterase, CVF/C3 (cobra venom factor/complement 3), 
fXaTx, fV aTx, hyaluronidase, L-amino acid oxidase and SVMP, displayed very little variation in their 
coding sequences, as did CRiSP & nerve growth factor (NGF). This is consistent with the mode of 
evolution generally adopted by large globular proteins [1]. In contrast, extensive variation was seen for 
3FTx, lectin, natriuretic, PLA 2 , kunitz and waprin toxin types. By examining the venom gland 
transcriptomes of a wide taxonomic range of neglected Australian elapid snake species, we have been 
able to gain deep insight into the molecular evolutionary history of major toxin classes. In addition to 
this, we have revealed that the toxic arsenals of the small Australian elapids, many of which are 
typically considered harmless to humans, are potentially similar in complexity to those of larger, 
"medically significant" species. 



Table 1. Diversity of toxin transcripts recovered from each elapid snake species. 



Species/Toxin 
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X 


X 


X 


X 


X 


X 
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X 






X 


X 






X 


X 


X 



fXaTx = factor Xa; fVaTx = factor Va; PLA2 = phospholipase A2; SVMP = snake venom metalloprotease; 



3FTx = three finger toxin. 
2.1. 3FTx 

3FTx are amongst the most abundant and well-studied components of elapid snake venoms [27]. 
The a-neurotoxic 3FTx from the venoms of Australian elapid snakes have been characterized into 
three groups: Types I, II (both also found in African and Asian elapids) and III (unique to Australian 
elapids [27]). Their cysteine arrangements and the number of residues present between cysteines [27] 
distinguish these three forms from one another. Type I (AKA short chain) a-neurotoxins are 
characterised by having lost the second and third cysteine residues present in the plesiotypic 3FTx 
form (leaving them with 8 cysteines), a change which may have resulted in a 100-fold increase in 
neurotoxicity [27]. Type II a-neurotoxins (AKA long chain) are characterised by having the same 
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eight cysteines, but with an additionally derived pair located between the fourth and fifth plesiotypic 
cysteine (third and fourth of the cysteines shared with Type I a-neurotoxins) [27]. The 
presence/absence of these two derived cysteines is key to the potency and specificity of the two 
a-neurotoxins. While both bind strongly to neuromuscular post-synaptic nicotinic acetylcholine 
receptors (nAChR), only the Type II can bind to neuronal nAChR [27]. 

Figure 2. Homology model depicting the locations of positively selected sites from various 
species, indicated by different colour codes. Multiple sequence alignment of Type II 
a-ntxs depicting the locations of positively selected sites is also presented. Representative 
sequences are from Brachyurophis roperi (1. GAHA01000012, 2. GAHA01000013, 
3. GAHAO 10000 16), Cacophis squamulosus (4. GAHB01000003, 5. GAHB01000008, 
6. GAHB0 1000008), Drysdalia coronoides (7. FJ752483, 8. FJ752485, 9. FJ752487), 
Hemiaspis signata (10. GAHF01000010, 11. GAHF0100001 1, 12. GAHF0 10000 14), 
Parasuta nigriceps (13. FJ790454, 14. FJ790448, 15. FJ790450), Vermicella annulata 
(16. GAHJ01000013, 17. GAHJ01000014, 18. GAHJ0 10000 15). Numerical IDs 
representing species lacking unique mutations are indicated by strikethrough. 
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The key functional sites for nAChR antagonizing activity in a-neurotoxins have been identified as 
being between residues 49 and 55 (between cysteines 3 and 4) in Type I a-neuro toxin and between 
residues 47 and 58 (between cysteines 3 and 6) in Type II a-neurotoxin [28]. The functional residues 
of the much smaller Type III a-neurotoxin remain to be elucidated. Variations at these functional 
sites, as well as variations in cysteine arrangements, are likely to have important consequences for the 
potency or affinity of these toxins and are thus of considerable interest in bioprospecting studies. 

156 distinct 3FTx sequences were recovered in this study (46 from P. modesta alone), making this 
by far the most diverse toxin type. Of particular interest amongst these are Type I and II a-neurotoxins 
with novel cysteine arrangements (Figures 2-6). A Type 1 isoform was recovered from 
C. squamulosus that possessed both the double cysteine (plesiotypic cysteine 1 and a novel cysteine) 
characteristic of and unique to some Australian elapid snake Type I a-neurotoxins, as well as an 
additional double cysteine (plesiotypic cysteine 2 paired with a novel cysteine). These toxins are part 
of an Australian elapid snake Type I a-neurotoxin clade with members that typically have an extra 
cysteine, so in this case the addition of a novel cysteine resulted in a unique ten cysteine arrangement, 
which may result in a novel folding pattern and activity. 

Figure 3. Molecular evolution of Type I (aka: short-chain) a-neurotoxins. 
Three-dimensional homology models of Type I a-neurotoxins from various species, 
depicting the locations of positively selected sites (Model 8, PP > 0.95, Bayes-Empirical 
Bayes approach) is presented here. Species are: (A) Brachyurophis roperi, (B) Cacophis 
squamulosus, (C) Drysdalia coronoides, (D) Hemiaspis signata, (E) Parasuta nigriceps 
and (F) Vermicella annulata. 
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Figure 4. Homology model depicting the locations of positively selected sites from various 
species, indicated by different colour codes. Multiple sequence alignment of Type II 
a-ntxs depicting the locations of positively selected sites is also presented. Representative 
sequences are from Acanthophis wellsi (1. GAGZ01000001, 2. GAGZ01000004, 
3. GAGZO 1000006), Brachyurophis roperi (4. GAHA01000003, 5. GAHA01000001, 
6. GAHA0 1000002), Drysdalia coronoides (7. FJ481928, 8. FJ752461, 9. FJ752459), 
Echiopsis curta (10. GAHD01000001, 11. GAHD01000004, 12. GAHDO 1000006), Furina 
ornata (13. GAHE01000001, 14. GAHE01000009, 15. GAHEO 10000 14), Hemiaspis 
signata (16. GAHF01000001, 17. GAHF01000005, 18. GAHFO 1000006), Suta fas data 
(19. GAHI01000001, 20. GAHIO 1000004), 21. Parasuta nigriceps FJ790457, Pseudonaja 
modesta (22. GAHH01000040, 23. GAHH01000045, 24. GAHH01000046, 
25. GAHH01000043, 26. GAHH01000042, 27. GAHH01000035) and Vermicella annulata 
(28. GAHJ01000009, 29. GAHJ01000010, 30. GAHJO 10000 16). Numerical IDs 
representing species lacking unique mutations are indicated by strikethrough. 
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Figure 5. Molecular evolution of Type II (aka: long-chain) a-neuro toxins. 
Three-dimensional homology models of Type II a-neurotoxins from various species, 
depicting the locations of positively selected sites (Model 8, PP > 0.95, Bayes-Empirical 
Bayes approach). Species are: (A) Acanthophis wellsi, (B) Brachyurophis roperi, 
(C) Drysdalia coronoides, (D) Echiopsis curta, (E) Furina ornata, (F) Hemiaspis signata, 
(G) Parasuta nigriceps and (H) Pseudonaja modesta. 
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Figure 6. Structural and functional evolution of Type III a-neuro toxins. Multiple sequence 
alignment of Type III a-ntxs depicting the locations of positively selected sites (Model 8, 
PP > 0.95, Bayes-Empirical Bayes approach) in various species of Australian elapids is 
presented here. Homology model depicting the locations of positively selected sites from 
various species, indicated by different colour codes, is also presented. Representative 
sequences are from Brachyurophis roperi (1. GAHAO 1000009, 2. GAHA0 10000 10, 
3. GAHAO 10000 11), Cacophis squamulosus (4. GAHB01000009, 5. GAHB01000010, 
6. GAHB0 10000 11), Furina ornata (7. GAHE01000022, 8. GAHE01000023, 
9. GAHE01000020, Pseudonaja modesta (10. GAHH01000009, 11. GAHH01000015, 
12. GAHH01000022) and Vermicella annulata (13. GAHJ01000001, 14. GAHJ01000003, 
15. GAHJ0 1000004). Numerical IDs representing species lacking unique mutations are 
indicated by strikethrough. 
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The Type II a-neurotoxins recovered included isoforms with newly evolved cysteines in addition to 
variants with cysteine deletions (Figures 4 and 5). An isoform recovered from B. roperi exhibiting 
deletion of the first of the newly evolved cysteines characteristic of Type II a-neurotoxins (located 
between plesiotypic cysteines 4 and 5 [27]), possessed an additional cysteine present near the 
C-terminus. This novel arrangement may have radical effects on disulphide bridging, thus exposing 
different residues upon the molecular surface of the toxin, with potential implications for relative 
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bioactivity or potency. In contrast, isoforms from P. modesta had deletions of both of the Type II 
a-neurotoxin characteristic cysteines. Sequences lacking these two cysteines have previously been 
reported from the venom of Pseudonaja textilis, Oxyuranus micro lepidotus and O. scutellatus [28]. 
Such a form from P. textilis has been shown in a prior study to be a blocker of both neuronal and 
neuromuscular nAChR [28]. This study did not, however, compare the neuronal binding affinity of this 
form to that of Type II a-neurotoxins with both these cysteines present. Thus change in relative 
potency or taxon specificity resulting from the loss of these cysteines remains unknown. The recovered 
diversity of molecular scaffolds is probably a result of the extreme influence of positive Darwinian 
selection experienced by a-neurotoxins [29]. 

3FTx sequence diversity extended beyond structural residues, with 3FTx recovered in this study 
differing from previously published 3FTx at key functional sites believed to be responsible for 
a-neurotoxic activity (Figures 2, 4 and 6). These include, amongst the Type I a-neurotoxin sequences, 
isoforms from B. roperi, C. squamulosus, F. ornata, H. signata and V. annulata. Amongst the Type II 
a-neurotoxin sequences, those with novel residues at functional sites include the aforementioned 
sequences (with deletion of one or both Type II characteristic cysteines) from B. roperi and 
P. modesta. While the key functional residues for the Type III a-neurotoxins have not been elucidated, 
the significant sequence diversity in the region of the loop-tips (Figures 6 and 7) suggests that 
significant variation in potency or specificity may exist for this toxin type as well. 

Figure 7. Molecular evolution of Type III a-neurotoxins. Three-dimensional homology 
models of Type III a-neurotoxins from various species, depicting the locations of 
positively selected sites (Model 8, PP > 0.95, Bayes-Empirical Bayes approach). Species 
are: (A) Furina ornata, (B) Pseudonaja modesta and (C) Vermicella annulata. 
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Phylogenetic analysis of 3FTx revealed that all sequences recovered in this study were placed into 
the Type I, II, or III a-neurotoxin clades previously characterised from Australian elapids [27]. Type 
III a-neurotoxins were recovered from P. modesta, which was expected as they have previously been 
characterised from other Pseudonaj 'a as well as Oxyuranus venoms [28,30,31]. Interestingly, they were 
also recovered from B. roperi, C. squamulosus, F. ornata and V. annulata, representing the first time 
this toxin type has been recovered from species outside the Pseudonaj ai 'Oxyuranus clade. Thus, it 
appears that this unique form was derived at the base of the Australian elapid snake radiation. 

As 3FTx sequences were by far the most numerous toxin sequences recovered, the molecular 
evolution of this toxin type in the venom system of Australian elapid snakes was further investigated. 
Integrative selection assessment using various methodologies (codeml site-specific models: M8, M2a, 
M3, MO; HyPhy: SLAC, FEL, REL, MEME, FUBAR, integrative analyses, branch-site REL and the 
evolutionary fingerprint analyses) revealed that most a-neurotoxins in the venoms of small Australian 
elapid snakes have evolved rapidly and episodically under the significant influence of positive 
selection (Figures 2-7; Supplementary Table 1). An exception is that of a V. annulata Type I 
a-neurotoxin, which was found to lack variation in coding sequences and appeared to have evolved 
under a regime of negative selection. Interestingly, Type II and Type III a-neurotoxins from this 
species accumulated variations rapidly under the influence of positive selection. Similarly, most Type I 
a-neurotoxin sequences recovered from H. signata were well conserved, while Type II a-neurotoxins 
from the same species were found to harbour large number of hypermutable sites. Since we only 
recovered three Type III a-neurotoxin sequences from B. roperi and C. squamulosus, selection 
assessment was not conducted for these sets. However, all Type III sequences in both these species 
shared a very high degree of sequence identity. In contrast, Type I and Type II a-neurotoxin 
sequences from B. roperi and Type I a-neurotoxin from C. squamulosus contained a greater amount 
of coding sequence variation. Thus, molecular evolutionary assessments revealed that the three kinds 
of a-neurotoxins in the venoms of small Australian elapid species have adopted differential 
evolutionary rates. Within each species, while one of the three forms of a-neurotoxin remains 
conserved, the other two accumulate significant variations, with the forms conserved or varied 
differing between species. 

We define focal mutagenesis as a phenomenon where the rapid accumulation of non-synonymous 
mutations under the influence of positive Darwinian selection in certain regions of the protein, such as 
the loops and the molecular surface of the toxin [29], has adaptive significance. Although, mutations 
occur randomly and the probability of a mutation occurring in structurally/functionally important and 
unimportant regions is theoretically equal, mutations in structurally/functionally important regions 
could result in the formation of destabilised and defective toxins. Since venom is energetically very 
expensive to produce [2], it would be a huge waste of resources for a venomous organism to secrete 
defective venom components. As a result, during the course of evolution, negative selection filters 
harmful mutations out of the populations and a greater majority of mutations are found in those regions 
that are structurally/functionally unimportant. Cysteine residues that allow the formation of disulphide 
bonds and in turn stabilise the toxin structure are one of such very important regions and hence are 
extremely well conserved in various venom-components. The mutation of surface chemistry under the 
influence of positive Darwinian selection, may not only enable toxins to non-specifically interact with 
novel prey receptors and generate a plethora of pharmacological reactions in the prey, but at the same 
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time might also allow them to escape/delay immunogenic response of prey animals. Not-surprisingly, 
focal mutagenesis has been documented in a diversity of venomous animal lineages [20,32-40]. 

Recently, we have demonstrated that various types of three-finger toxins in elapid venoms have 
adopted RAVER, and accumulate hypermutable sites in their loops and on the molecular surface, 
while conserving structurally and functionally important residues [29]. Mapping of hypermutable sites 
on the alignments and three-dimensional homology models of various 3FTx functional forms 
recovered from various small Australian elapids, revealed a number of interesting evolutionary 
phenomena: (i) a-neurotoxins from various species have accumulated variations in similar regions of 
the toxin, highlighting the role of focal mutagenesis in the evolution of these genes; (ii) certain species 
have harboured unique hypermutable sites, which might account for the differences in efficacy and/or 
potency of these toxins; (iii) residues in a-neurotoxins that have been elucidated as structurally or 
functionally important in other species were found to be extremely well conserved in most Australian 
elapids, while most hypermutable sites were found in regions not known to be important for function 
or structural stability; (iv) most hypermutable sites were concentrated in the loops of the toxins and (v) 
sequences from different species have different sets of residues that exhibit conservation, indicating 
that the diet of these snakes could have influenced the evolution of a-neurotoxin gene. Thus, 
molecular evolution assessment not only revealed that a-neurotoxins in the venoms of small 
Australian elapid snakes have adopted differential rates of evolution and focal mutagenesis, but also 
suggested that there may be a correlation between the feeding ecology of a species and the evolution of 
its a-neurotoxin-encoding gene. The diversity of residues found in functionally important sites of 
a-neurotoxins from species with differing feeding ecologies indicates that the evolution of these toxins 
in small Australian elapid snakes may have been driven by dietary specialisation. Investigations of the 
venom proteome of these snakes may be illuminating in this regard as may pharmacological 
investigation of the potential prey-specificity of individual a-neurotoxin isoforms and/or clades. 

2.2. Lectin 

The lectin sequences recovered in this study considerably expand the number of lectin sequences 
recovered from the venoms of elapid snakes. The plesiotypic form of lectins from reptile venom has 
been demonstrated to possess a "carbohydrate specificity triad" of functional residues that confers 
either galactose-binding (QPD motif — [41]) or mannose-binding (EPN motif — [42]) mediated 
haemagglutination activity, with EPN as the plesiotypic state. The majority of lectins recovered from 
this study had the EPN motive. However, sequences with the galactose-binding QPD motif were 
recovered from Cacophis squamulosus and Pseudonaja modesta. In addition to these two 
functionally-characterised motifs, sequences were recovered in this study that had the novel motifs 
KPG {Acanthophis wellsi) or YRH (Cacophis squamulosus) at the carbohydrate binding site (data not 
shown). These divergent sequences may prove to have novel activities. 

2.3. Natriuretic Peptides 

The venom gland transcriptomes of several species of elapid snake examined in this study 
(e.g., E. curtus, H. signata) contained numerous distinct isoforms of natriuretic peptide, suggesting that 
this toxin type may be a more important component of the venom of some species than hitherto 
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realised. 

C-type natriuretic peptides (CNP) from snake venom are potent vasodilators and hypotensive agents 
(c.f. [43,44]). In addition, the precursor contains multiple proline rich peptides, one of which was 
utilised in the development of the multibillion-dollar drug Captopril [45]. For this reason they are 
considered to be amongst the most promising of all animal toxins for pharmaceutical 
bioprospecting [4-7,46]. Despite this high level of interest, very few complete sequence precursors of 
snake venom natriuretic peptide precursors have been recovered. The final processed peptide of the 
natriuretic domain has been sequenced from Austrelaps superbus; Cryptophis nigrescens; 
Hoplopcephalus stephensii; Notechis scutatus; Oxyuranus micro lepidotus; Oxyuranus scutellatus; 
Pseudechis australis; Pseudechis porphyriacus; Pseudonaja textilis; and Tropidechis 
carinatus [43,47,48] while only a single complete precursor had been sequenced prior to this study 
(from O. scutellatus (Uniprot accession P83228)). Due to this low level of sequencing, our 
understanding of the evolutionary history of this toxin type was far from complete. 

In this study, in addition to retrieving natriuretic sequences, which were similar to the previously 
characterised precursors from the elapid snake Oxyuranus scutellatus (Uniprot accession P83228), 
isoforms were recovered that lacked the derived C-terminal tail typical of elapid natriuretic toxin 
variants (Figure 8). Instead, these toxins resembled viperid and "non-front-fanged" (NFF) advanced 
snake natriuretic homologues in lacking the tail and possessing glycine rich regions characteristic of 
the aforementioned forms. Phylogenetic analyses demonstrated the relationship of these novel elapid 
snake sequences to the tail-less natriuretic sequences recovered from viperid and colubrid snakes 
(Figure 9). 

Natriuretic peptides in elapid snake venoms have previously been inferred to have been derived 
from ANP (atrial natriuretic peptide) or BNP (brain natriuretic peptide) due to the presence of the 
C-terminal tail, while those of viperid snakes were said to be derived from CNP (c-type natriuretic 
peptide) due to lack of the tail [1,44,49-51]. However, recent studies have shown that all snake 
natriuretic toxins form a monophyletic group within the CNP clade [15,22]. Despite this, the 
evolutionary linkage between elapid and viperid sequences has remained enigmatic. It has been 
suggested that presence of the C-terminal extension is an apotypic (derived) condition of the peptide 
and that the tail-less form, previously only recovered from viperid snakes, is the plesiotypic 
condition [22]. The recovery in this study, for the first time, of tail-less forms from elapid snakes 
(C. squamulosus and S. fasciata) supports this hypothesis, particularly as these sequences also possess 
the glycine-rich region of the propeptide that is characteristic of the tail-less viperid snake forms. In 
phylogenetic analyses performed in the present study, the tail-less elapid snake forms grouped with 
tail-less forms from viperid and colubrid snakes and the sole published viperid snake precursor (from 
C. cerastes) possessing the C-terminal extension grouped with elapid and lamprophiid snake 
sequences that share the tail (Figure 9). These results add further strength to the hypothesis that the 
C-terminal extension represents a derivation of this toxin type and thus the action of the tailed forms 
on GC-A [43] is not indicative of descent from the atrial natriuretic peptide (ANP) but instead 
represents a remarkable case of convergence for receptor-specific targeting [22]. The C-terminal 
extension appears to have evolved in elapid snake venom natriuretic peptides to confer greater affinity 
for GC-A (guanylate cyclase A) and NPR-C (natriuretic peptide receptor C) receptors [43] and greater 
resistance to proteolytic cleavage from the receptors [52]. Recently a study was published in which a 



Toxins 2013, 5 



2634 



natriuretic peptide possessing the C-terminal extension was sequenced from the venom of the 
rattlesnake Crotalus oreganus abyssus [53]. The published sequence represents only the final 
processed natriuretic peptide amino acid sequence and does not provide precursor information. As all 
snake sequences group into a monophyletic clade nested within the C-type natriuretic peptide (CNP) 
clade of non-venom precursors (Figure 9) [15,22], the venom form is clearly of CNP ancestry and the 
C-terminal tail evolved soon after recruitment of this form to the venom system. 

Figure 8. Sequence alignment of natriuretic peptides. (1). P68515 Bothrops insularis, 
(2). K4J3K2 Azemmiops feae, (3). K4IT20 Azemmiops feae, (4). A8YPR6 Echis ocellatus, 
(5). Q09GK2 Philodryas olfersii, (6). GAHIO 10000 13 Suta fasciata, (7). P83228 
Oxyuranus scutellatus, (8). GAHIO 10000 16 Suta fasciata, (9). A8YPR9 Cerastes cerastes 
Post-translationally cleaved peptides in shaded in gray. 
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Figure 9. Phylogenetic reconstruction of the molecular evolutionary history of natriuretic 
peptides. Non-toxin outgroup sequences (P23582 and P55207) not shown. Representative 
sequences obtained in this study are shown in red. Node labels indicate 
posterior probabilities. 
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Recently, a proline-rich peptide was isolated from the venom of Dendroaspis angusticeps [54] 
which was described as being of an unknown peptide or protein type. The results of this study reveal 
that this proline-rich peptide is part of the propeptide region of the natriuretic peptide precursor, a fact 
missed by the aforementioned Dendroaspis study despite an Oxyuranus precursor being available that 
showed significant similarity to this peptide in an early stretch of the propeptide region. Such a region 
was preserved in the sequences obtained in this study. This region is posttranslationally cleaved from 
the precursor and undergoes dynamic evolution facilitating neofunctionalisation. It is where the 
bradykinin-potentiating peptides (of Captopril fame) evolved, and is also the source of the novel 
neurotoxic peptides from Azemiops feae and Psammophis mossambicus [32,55]. Thus, it represents an 
unexplored source of novel peptides with potentially useful activities for drug design and development. 

2.4. PLA 2 

PLA 2 toxins are an important component of the venom of many Australian elapid snakes and 
typically exhibit presynaptic neurotoxic activity, myotoxic activity, or both, although some forms 
exhibit antiplatelet activity (c.f. [14]). Almost all previous sequences have come from the well-studied 
larger species, despite toxic effects of venoms from the smaller snakes indicating they are rich in PLA 2 
toxins. For example Furina (Glyphodon) tristis venom has been demonstrated to be presynaptic ally 
and myotoxically active [56], indicating the presence of PLA 2 toxins, a result consistent with 
components detected in the venom of this species by LC/MS (liquid chromatography mass 
spectrometry) that had PLA 2 characteristic retention times and molecular weights [57]. In addition, the 
venoms of two species of small elapid snake from the genus Suta have been pharmacologically 
characterised as potently neurotoxic (S. suta) or myotoxic (S. punctata) [58]; mass spectrometry of 
these venoms detected components consistent with PLA 2 as well as three-finger toxins [57]. The 
Kuruppu 2007 study could not rule out the presynaptic mode of action for the observed neurotoxicity, 
while the myotoxicity was almost certainly caused by the action of PLA 2 . These results demonstrate 
that the venoms of small Australian elapid snakes contain similar toxins to those of their larger 
cousins, and are potentially capable of delivering serious bites. This potential has been made 
alarmingly clear in the case of Furina tristis and Suta suta, both of which have been responsible for 
life-threatening envenomations in which the neurotoxic action of the respective venoms was poorly 
reversed by polyvalent antivenom ([10], BG Fry, personal observations). The results of the present 
study considerably augment the number of sequenced Type II PLA 2 full length sequences from elapid 
snakes. The recovery of PLA 2 sequences from 8 of the 1 1 libraries confirms that PLA 2 is one of the 
most important and highly expressed toxin types in the venom of Australian elapid snakes. The new 
sequences recovered in the present study also provide further insight into the molecular evolutionary 
history of this clinically significant toxin type (Figure 11). 

Several previously published PLA 2 sequences are known to be components of multimeric 
neurotoxin complexes (e.g., taipoxin gamma chain [59]). These sequences feature a novel cysteine 
arrangement with the insertion of two additional cysteines between plesiotypic cysteines 1 and 2. 
Multimeric presynaptic PLA 2 have long been known from the venoms of Oxyuranus sp. and 
Pseudonaja sp and are amongst the most potent neurotoxins known. Such complexes have recently 
been reported to be present in the venoms of Acanthophis sp. [60-63]. 



Toxins 2013, 5 



2637 



In the present study, sequences homologous with all three chains of taipoxin were recovered from 
A. wellsi and such transcripts were also present in the S. fasciata library (Figure 10). This supports a 
previous hypothesis that they were likely widespread in the venoms of Australian elapids [64]. 
Phylogenetic analyses showed that these Acanthophis and Suta sequences grouped into clades with 
each of the respective taipoxin subunits (a, P and y) from Oxyuranus scutellatus (Figure 11). 
Intriguingly, the P. modesta library did not include any of the subunits of the potent presynaptic 
neurotoxic PLA 2 complex 'textilotoxin' present in the venom of the congeneric Pseudonaja textilis 
(Figure 11). 

Figure 10. Sequence alignment of 'taipoxin/paradoxin'-like presynaptic complex subunits: 
a-subunit (1). Q45Z43 Oxyuranus micro lepidotus, (2). Q45Z48 Oxyuranus scutellatus, 
(3). GAGZO 1000028 Acanthophis wellsi, (4). A6MFM9 Rhinoplocephalus nigrescens, 
(5). GAHI01000025 Suta fasciata, (6). B5G6G1 Tropidechis carinatus, P-subunit (7). 
Q45Z46 Oxyuranus micro lepidotus, (8). Q45Z53 Oxyuranus scutellatus, 
(9). GAGZO 1000024 Acanthophis wellsi, (10). GAHI0 1000027 Suta fasciata and 
y-subunit (11). Q4VRI6 Oxyuranus scutellatus, (12). GAGZ01000027 Acanthophis wellsi, 
(13). Q9PUG7 Austrelaps superbus, (14). GAHI01000030 Suta fasciata. 
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Figure 11. Phylogenetic reconstruction of the molecular evolutionary history of snake 
venom Type I phospholipase A 2 toxins. Non-toxin outgroup sequences (Q8JFB2 and 
Q8JFG2) not shown. Representatives of sequences obtained in this study are shown in red. 
Node labels indicate posterior probabilities. 
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2.5. Kunitz and Waprin 

An interesting picture emerges from the results of this study in relation to the molecular 
evolutionary histories of the kunitz and waprin peptides. In addition to the single domain kunitz and 
waprin sequences, dual domain (kunitz -kunitz) and fused dual domain (kunitz-waprin) precursors were 
also recovered (the latter representing the first time fused kunitz-waprin precursors have been 
recovered from elapid snakes) (Figure 12). It is noteworthy that a mono-domain waprin sequence 
recovered from Denisonia devisi was highly similar to a previously published sequence from the 
natricid NFF advanced snake Rhabdophis tigrinus [19], and shared an identical signal peptide. Dual 
domain (kunitz -kunitz) kunitz peptide sequences were recovered from the venom gland transcriptomes 
of Furina ornata (fragment), Hoplocephalus bungaroides and Hemiaspis signata and fused dual 
domain kunitz-waprin toxin sequences were recovered from Cacophis squamulosus and Suta fasciata. 
Signal peptides differed considerably between the kunitz-kunitz dual domain precursors and the single 
domain kunitz precursors, the single domain waprin peptides and kunitz-waprin hybrid precursors (all 
of which share the same signal peptide in elapid snake forms [65] other than the novel waprin 
sequences fromZ). devisi mentioned above-Figure 12). 

The fact that both dual-domain precursors of these toxin types are widespread and that the 
mammalian physiological homologues also have dual-domain precursors for both peptide types [66] 
strongly suggests that the plesiotypic form of each toxin type is a dual-domain precursor, consistent 
with phylogenetic analyses of the relative position of the dual-domain form [20,33]. Thus the 
mono-domain forms are secondary derivations resulting from the loss of a domain. In both toxin types 
it appears to be the first domain that has been deleted. It is worth noting, pending proteomic 
investigation, that it remains unknown whether or not the kunitz-waprin fused toxin is 
posttranslationally processed into two separate toxins or if it acts as a complex of hitherto 
uncharacterised activity. 

The origin of the shared signal peptide of mono-domain kunitz and waprin peptides, as well as 
fused kunitz-waprin forms, is enigmatic. It has been suggested that kunitz and waprin arose from 
duplication of the same ancestral gene and that subsequent diversification occurred only in the 
toxin-encoding region, thus preserving the sequence homology of the signal peptides [67,68]. This 
hypothesis is in conflict with the fact that the ancestral dual-domain encoding toxin forms of each of 
these peptides types possess distinct signal peptides. Mono-domain waprin peptides sequenced from 
colubrid snakes also possess a different signal peptide to previously published mono-domain forms 
from elapid snakes, although in the present study a sequence with high degree of similarity (including 
the same signal peptide) to a sequence from the natricid snake Rhabdophis tigrinus [19] was recovered 
from the elapid snake Denisonia devisi. Only two (mono-domain) kunitz precursors have been 
recovered from colubrid snakes to date, one of which (Telescopus dhara) possesses the shared signal 
peptide, whilst the other (Philodryas olfersii) possesses yet another distinct signal peptide [19]. 
Mono-domain kunitz from viperid snakes possess the shared signal peptide (c.f. [69]). As yet, no 
mono-domain waprin peptide precursors have been sequenced from viperid snakes. It is apparent that a 
gene-splicing event took place at some point in the molecular evolutionary history of these toxins but 
the origin of the shared derived signal peptide is unclear. Phylogenetic analysis of kunitz sequences 
recovered in this study, along with previously published sequences, revealed that the dual 
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kunitz-domain encoding precursors group together as a basally-divergent clade, further strengthening 
the hypothesis that the dual-domain precursor is the plesiotypic form of this toxin type. 

Figure 12. Sequence alignment of precursors encoding: dual-domain kunitz (1). B2BS84 
Austrelaps labialis, (2). GAHGO 1000009 Hoplocephalus bungaroides; mono-domain 
kunitz (3). GAGZ0 10000 19 Acanthophis wellsi, (4). GAGZ0 10000 17 Acanthophis wellsi, 
(5). GAHB0 10000 16 Cacophis squamulosus, (6). GAHD0 10000 11 Echiopsis curta, 
(7). GAHGO 1000008 Hoplocephalus bungaroides, (8). GAHH0 1000051 Pseudonaja 
modesta, (9). GAHI0 10000 10 Suta fasciata; dual-domain waprin (10). A7X4K1 
Philodryas olfersii; mono-domain waprin (11). GAHC0 1000021 Denisonia devisi, 
(12). A7X4J4 Rhabodophis tigrinus, (13). A7X4K7 Philodryas olfersii, (14). A7X4I7 
Thrasops jacksonii, (15). B5G6H4 Notechis scutatus, (16). B5G6G8 Oxyuranus 
scutellatus; kunitz-waprin fusion (17). D3U2B9 Sistrurus catenatus edwardsii, 
(18). D3U0D3 Sistrurus catenatus tergeminus, (19). GAHB01000034 Cacophis 
squamulosus, (20). GAHI0 1000009 Suta fasciata. 
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The function of waprin in snake venom is almost entirely unknown. Only one component, from the 
venom of O. microlepidotus, has been tested and this toxin, named "omwaprin", was found to exhibit 
weak antimicrobial activity [70]. However, phylogenetic analysis of sequences recovered in this study 
revealed an interesting level of sequence diversity. Although the majority of elapid sequences grouped 
together, a sequence recovered from C. squamulosus in this study was found to be basally-divergent 
from the main clade of elapid sequences while sequences recovered from A. wellsi, D. devisi and V. 
annulata displayed affinity to previously published sequences from colubrid snakes and African elapid 
snakes. The clustering of the fused kunitz-waprin precursors within the main clades of waprin or 
kunitz toxin sequences in the respective analyses suggests that the mono-domain toxin forms of these 
peptides had already been derived prior to their fusion. Nonetheless, the presence of the fused toxin in 
both viperid and elapid snakes, which are the most divergent pair of colubroid snake families, suggests 
that the kunitz-waprin splicing is relatively ancient, and occurred at the base of the 
Colubroidea radiation. 

2.6. Procoagulant Toxins and Clinical Implications 

As well as theoretical implications regarding molecular evolution, our discoveries may prove 
important in considerations of the potential clinical effects of bites from these species, such as our 
sequencing of fXaTx from H. signata and H. bungaroides, which is consistent with coagulopathic 
effects produced by the venom of these snakes [11,71,72]. Of interest is the inability in the current 
study to recover fXaTx sequences from the venom gland transcriptome of P. modesta; this is 
significant since high expression levels of fXaTx and fVaTx are characteristic of the venom of other 
members of the genus Pseudonaja [13,14]. Although it is possible that fXaTx is expressed in very low 
levels and that the sampling in this study simply failed to detect it, it is also possible, as previously 
suggested [73] that the venom of this species, unusually for a member of its genus, lacks any 
procoagulant component to its venom. The recovery of fVaTx from P. modesta in the present study, 
however, raises further interesting questions. A rate-limiting step in the action of fXaTx is that it must 
form a 1:1 complex with activated factor V in the bitten animal's blood [74]. fVaTx is present in the 
venoms of species Pseudonaja and all studied species of Oxyuranus and, as it has not been detected in 
the venoms of any other species, it is believed to have been recruited into the venom system of the 
common ancestor of both genera, which form a monophyletic clade [14,75]. Not only is this toxin 
responsible for the potentially fatal venom-induced consumption coagulopathy (ViCC) often suffered 
by victims of bites from these snakes [76], it has also been hypothesised that this toxin played an 
important role in the evolutionary transition of snakes in this clade from a diet of reptiles and frogs 
(common to most Australian elapid snakes including the species included in this study) to a diet 
specialising in mammalian prey [19,77]. In the venom of these snakes, fVaTx forms a complex with 
fXaTx creating a "complete prothrombin activator" [78]. The function of the venom form of fVaTx is 
to act as a cofactor for fXaTx that doesn't require proteolytic cleavage to function and is resistant to 
the regulatory activity of the anticoagulant protein C [79]. The venom form of fVaTx is therefore a 
cofactor for fXaTx that is faster acting and harder to stop than its endogenous counterpart. Why 
P. modesta might express this cofactor in the absence of fXaTx is a mystery. The absence of fXaTx 
may, however, relate to the diet of this species — P. modesta, uniquely for a member of its genus, 
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appears to prey exclusively upon reptiles [80]. It is also worth noting that only very low expression 
levels of fVaTx were detected in the P. modesta library, which indicates that its expression is probably 
down-regulated in this species relative to other Pseudonaja sp. In order to absolutely confirm the 
absence of fXaTx in the venom system of P. modesta and in order to shed light on the recruitment date 
of this important toxin, it will be necessary in a future study to use PCR amplification with primers to 
facilitate the detection of sequences which may be expressed at very low levels. The radically different 
venom profile of P. modesta however is indicative of a potential failure of CSL Brown Snake 
antivenom to neutralise envenomation effects as the primary antibodies in the antivenom are selected 
towards the most dominant toxins in the immunising venom (Pseudonaja textilis), and thus would 
likely be directed predominantly against the fXaTx-fVaTx complex and not against 3FTx which is the 
dominant toxin type in P. modesta venom. This is crucially important as 3FTx been hypothesised to 
poorly immunogenic despite the high toxicity [81]. 

2. 7. Influence of Prey Preference on Venom Gland Transcriptome 

There is abundant evidence that prey preference drives the evolution of venom composition and 
toxicity in snakes [4,20,37,77,82-84]. It is perhaps surprising, therefore, that a reptile specialist like S. 
fasciata should have such a complex venom gland transcriptome, but it remains to be demonstrated 
whether or not the venom proteome is similarly complex or if it is complex but dominated by one 
particular toxin type or particular isoforms within a toxin type. Even the venom gland transcriptome of 
V. annulata, a species with venom that has been previously been inferred (via mass spectrometry 
analysis) to be dominated by a single toxin class (3FTx — [11]), contained a diversity of other toxin 
types. V. annulata is a burrowing species with a highly specialised diet; apparently feeding only on 
blind snakes [85]. One fairly intuitive hypothesis is that species with highly specialised diets are likely 
to have simpler venom compositions than those that feed on a wide range of prey. The results from the 
aforementioned study [11] supported this hypothesis, however the results of the present study, in 
which precursors for 12 distinct toxin types were recovered from the transcriptome of V. annulata, 
appear to contradict it. One possible explanation is that V. annulata is descended from a species with a 
broader diet and hence a more diverse venom composition. If that is the case, it is likely that 
microRNA is blocking the translation of these precursors into toxins [86], as the snake, with its derived 
diet of blind snakes, no longer requires them. Further proteomic investigation of the venom of this 
species will confirm whether or not its venom is indeed streamlined and subsequent pharmacological 
investigation of the dominant toxins will confirm whether or not they exhibit taxon-specific potency 
variation. Confirming these properties of the venom and toxins would strongly support the hypothesis 
that V. annulata and its prey are engaged in a co-evolutionary predator-prey chemical arms race. It 
should be noted that a complex toxin cocktail might also be helpful for defensive purposes. 

2.8. Brachyurophis roperi: A Unique, Oophagous Burrowing Elapid Snake 

The preponderance of unique sequences within the B. roperi venom gland transcriptome is 
something of a mystery. B. roperi are believed to be specialist egg-eaters and have the most specialised 
dentition for oophagy of any terrestrial Australian elapid snake (Scanlon and Shine, 1988). In 
oophagous marine elapid snakes, a degeneration of the venom system has been noted, along with 
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deletions of functional residues in 3FTx and accumulation of deleterious mutations in PLA 2 [87,88]. It 
is possible that a similar scenario is taking place within the venom system of B. roperi. Unlike the 
oophagous sea snakes Aipysurus eydouxii, A. mosaicus and Emydocephalus annulatus, however, B. 
roperi appears to retain fangs that are similar in size to those of non-specialist congeners that feed on 
both eggs and lizards [89]. In addition, a bite from B. roperi resulted in localised pain and 
inflammation, which lasted for several hours, suggesting that this species is still in possession of a 
functional venom system able to cause muscle pain radiating up an entire arm lasting 12 hours (BG 
Fry, personal observations). 

Interestingly, Type I a-neuro toxins, which were found to be extremely well conserved in most 
Australian elapids, accumulated a large number of hypermutable sites in B. roperi; particularly in 
regions that harbour residues that confer a-neurotoxins in other species the ability to target muscular 
al -nicotinic acetylcholine receptors (nAChRs). Similarly, Type II a-neurotoxins in this species 
appeared to have poorly conserved functionally important residues, which confer Type II 
a-neurotoxins in other species the ability to target muscular al and neuronal al nAChRs. Thus, B. 
roperi Type I and II a-neurotoxins have likely lost their ability to target nAChRs of the prey. Hence, 
the observed pathogenesis could have been the result of other venom-components, including the 
relatively well-conserved Type III a-neurotoxins, or non-specific immunogenic reactions. In contrast, 
a-neurotoxin genes in other Australian elapids tend to accumulate variations largely in structurally and 
functionally unimportant regions. Hence, these regions are likely involved in mounting immunogenic 
reactions in the prey. 

3. Experimental Section 

3.1. Study Species 

The widest possible phylogenetic and ecological diversity amongst the smaller forms of Australian 
elapid snake were chosen for this study. Unlike the larger, more commonly studied species of 
Australian elapid snakes, which typically have generalized diets, the species chosen for this study feed 
predominantly on other reptiles (or their eggs, in the case of Brachyurophis roperi) except for 
Denisonia devisi, which feeds almost exclusively on frogs [12,80,85,90-95]. As well as reptile and 
frog specialists, this study also examines species with divergent ecologies including: burrowing 
species (Brachyurophis roperi and Vermicella annulata); a semi-arboreal species (Hoplocephalus 
bungaroides); ambush hunters (Acanthophis wellsi and Denisonia devisi); nocturnal foragers 
(Cacophis squamulosus and Furina ornata); arid zone inhabitants (Acanthophis wellsi and Suta 
fasciata); and marsh dwellers (Echiopsis curta and Hemiaspis signata). Also examined were 
anomalous and unstudied members of otherwise well studied genera: Pseudonaja modesta (the 
smallest Pseudonaja, a reptile specialist in a genus of generalists) and Acanthophis wellsi (an arid zone 
death adder). To date, no toxin sequences have previously been retrieved/published from any of the 
species examined in the present study. 

Species and collection localities: Acanthophis wellsi — Millstream, Western Australia; 
Brachyurophis roperi — Kununurra, Western Australia; Cacophis squamulosus — Mt Glorious, 
Queensland; Denisonia devisi — Glenmorgan, Queensland; Echiopsis curta — Perth, Western Australia; 
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Furina ornata — Kununurra, Western Australia; Hemiaspis signata — Mt Glorious, Queensland; 
Hoplocephalus bungaroides — Sydney, New South Wales; Pseudonaja modesta — Sandstone, Western 
Australia; Suta fasciata — Sandstone, Western Australia. Vermicella annulata — Mt. Glorious, 
Queensland. All specimens were adult males. 

3.2. Trans criptome Sequencing 

Total RNA was extracted from venom glands using the standard TRIzol Plus method (Invitrogen). 
Extracts were enriched for mRNA using standard RNeasy mRNA mini kit (Qiagen) protocol. mRNA 
was reverse transcribed, fragmented and ligated to a unique 10-base multiplex identifier (MID) tag 
prepared using standard protocols and applied to one PicoTitrePlate (PTP) for simultaneous 
amplification and sequencing on a Roche 454 GS FLX+ Titanium platform (Australian Genome 
Research Facility). As each plate contained mRNA samples from multiple species, automated 
grouping and analysis of sample-specific MID reads informatically separated sequences from the other 
transcriptomes on the plates, which were then post-processed to remove low quality sequences before 
de novo assembly into contiguous sequences (contigs) using v 3.4.0.1 of the MIRA software program. 

Full assembly parameters are available in Supplementary File 1 , assembly details in Supplementary 
Tables 2 and 3 and nucleotide sequences are available in Supplementary File 2 as well from Genbank: 
Acanthophis wellsi (GAGZ01 00000 1-GAGZ0 1000032), Brachyurophis roperi (GAHA0 100000 1 - 
GAHA0 1000031), Cacophis squamulosus (GAHB01 00000 1-GAHB0 1000034), Denisonia devisi 
(GAHC01 00000 1-GAHC0 1000021), Echiopsis curta (GAHD01 00000 1-GAHD0 1000032), Furina 
ornate (GAHE01 00000 1-GAHE0 100003 8), Hemiaspis signata (GAHF01 00000 1-GAHF0 1000041), 
Hoplocephalus bungaroides (GAHG01 00000 1-GAHG0 1000046), Pseudonaja modesta 
(GAHH01 00000 1-GAHH0 1000073), Suta fasciata (GAHI01 00000 1-GAHI0 1000020), and Vermicella 
annulata (GAHJ01 00000 1-GAHJ0 1000024). All raw reads have been deposited in the NCBI 
Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra/) with the accession numbers of: 
SRR768900 Acanthophis wellsi; SRR768902 Brachyurophis roperi; SRR768909 Cacophis 
squamulosus; SRR768910 Denisonia devisii; SRR768911 Echiopsis curta; SRR768912 Furina 
ornata; SRR768914 Hoplocephalus bungaroides; SRR768915 Pseudonaja modesta; SRR768916 Suta 
fasciata; SRR768917 Vermicella annulata. Assembled contigs were processed using CLC Main Work 
Bench (CLC-Bio) and Blast2GO bioinformatic suite to provide Gene Ontology, BLAST and 
domain/Interpro annotation. The above analyses assisted in the rationalisation of the large numbers of 
assembled contigs into phylogenetic 'groups' for detailed phylogenetic analyses outlined below. 

3.3. Phylogenetics 

Phylogenetic analyses were performed to allow reconstruction of the molecular evolutionary history 
of each toxin type for which transcripts were bioinformatically recovered. Toxin sequences were 
identified by comparison of the translated DNA sequences with previously characterised toxins using a 
BLAST search of the UniProtKB protein database. Molecular phylogenetic analyses of toxin 
transcripts were conducted using the translated amino acid sequences. Comparative sequences from 
other venomous reptiles and physiological gene homologs identified from non-venom gland 
transcriptomes were included in each dataset as outgroup sequences. To minimize confusion, all 
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sequences obtained in this study are referred to by their Genbank accession numbers 
(http://www.ncbi.nlm.nih.gov/sites/entrez?db=Nucleotide) and sequences from previous studies are 
referred to by their UniProtKB accession numbers (www.uniprot.org). Resultant sequence sets were 
aligned using CLC Mainbench. When presented as sequence alignments, the leader sequence is shown 
in lowercase and cysteines are highlighted in black. > and < indicate incomplete N/5' or C/3' ends, 
respectively and * used to indicate the end of a sequence. Datasets were analysed using Bayesian 
inference implemented on MrBayes, version 3.2.1 using lset rates=invgamma with prset 
aamodelpr=mixed, which enables the program to optimize between nine different amino acid 
substitution matrices. The analysis was performed by running a minimum of 1 x 10 7 generations in 
four chains, and saving every 100th tree. The log-likelihood score of each saved tree was plotted 
against the number of generations to establish the point at which the log likelihood scores reached their 
asymptote, and the posterior probabilities for clades established by constructing a majority-rule 
consensus tree for all trees generated after completion of the burn-in phase. 

3.4. Selection Analyses 

We evaluated the influence of natural selection on various types of three-finger toxins using 
maximum-likelihood models [96,97] implemented in CODEML of the PAML package of 
programs [98]. We employed site-specific models that estimate positive selection statistically as a 
non-synonymous-to-synonymous nucleotide-substitution rate ratio (co) significantly greater than 1 . We 
compared likelihood values for three pairs of models with different assumed co distributions as no 
a priori expectation exists for the same: MO (constant co rates across all sites) versus M3 (allows co to 
vary across sites within V discrete categories, n > 3); Mia (a model of neutral evolution) where all 
sites are assumed to be either under negative (co < 1) or neutral selection (co = 1) versus M2a (a model 
of positive selection) which in addition to the site classes mentioned for Mia, assumes a third category 
of sites; sites with co > 1 (positive selection) and M7 (Beta) versus M8 (Beta and co), and models that 
mirror the evolutionary constraints of Ml and M2 but assume that co values are drawn from a beta 
distribution [99]. Only if the alternative models (M3, M2a and M8: allow sites with co > 1) show a 
better fit in Likelihood Ratio Test (LRT) relative to their null models (MO, Mia and M7: do not allow 
sites co > 1), are their results considered significant. LRT is estimated as twice the difference in 
maximum likelihood values between nested models and compared with the % distribution with the 
appropriate degree of freedom-the difference in the number of parameters between the two models. 
The Bayes empirical Bayes (BEB) approach [100] was used to identify amino acids under positive 
selection by calculating the posterior probabilities that a particular amino acid belongs to a given 
selection class (neutral, conserved or highly variable). Sites with greater posterior probability 
(PP > 95%) of belonging to the 'co > 1 class' were inferred to be positively selected. 

Single Likelihood Ancestor Counting (SLAC), Fixed-Effects Likelihood (FEL) and Random 
Effects Likelihood (REL) implemented in HyPhy [101] were employed to detect sites evolving under 
the influence of positive and negative selection. Fast, Unconstrained Bayesian AppRoximation 
(FUBAR) [102] was also employed to detect sites evolving under the influence of pervasive 
diversifying and purifying selection. Mixed Effects Model Evolution (MEME) [103] was used to 
efficiently detect episodically diversifying sites. To clearly depict the proportion of sites under 
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different regimes of selection, an evolutionary fingerprint analysis was carried out using the 
evolutionary selection distance (ESD) algorithm [104] implemented in datamonkey Branch-sire REL 
tests [105] were employed to identify lineages evolving under episodic bursts of adaptive selection 
pressures. For all the aforementioned selection analyses, phylogenetic trees were built using the 
maximum-likelihood methodology implemented in PhyML 3.0 [106]. Node support was evaluated 
with 1,000 bootstrapping replicates. 

3.5. Structural Analyses 

To depict the natural selection pressures influencing the evolution of various three-finger toxins, we 
mapped the sites under positive selection on the homology models created using Phyre 2 
webserver [107]. Pymol 1.3 [108] was used to visualize and generate the images of homology models. 
Consurf webserver [109] was used for mapping the evolutionary selection pressures on the 
three-dimensional homology models. 

4. Conclusion 

Whilst the findings of the present study have significantly increased our knowledge of the 
molecular evolution of Australian elapid venom, little can be ascertained about the specific roles of 
these toxins in the lifestyles of these snakes without further work being conducted. Thus, future work 
must be undertaken to characterise the venom proteomes of the 1 1 species examined in this study. As 
well as shedding light on methodological questions (e.g., whether or not a complex transcriptome 
always equals a complex proteome), proteomic characterisation of these venoms will allow us to 
further investigate correlations between diet and venom composition (c.f. [82]), gauge their relative 
toxicity in different animal models, and more accurately assess the potential danger of these snakes to 
human bite victims. Alongside proteomic characterisation, it is important that the extent and basis of 
any cross-reactivity of these venoms with currently available antivenoms and the snake venom 
detection kit (SVDK) [1 10] be ascertained. 

Additional questions regarding the molecular evolution of these toxins remain to be answered, 
including the timing of the recruitment event of the fXaTx into the venom arsenal of Australian 
elapids. This question may be answered through the use of primer-driven PCR amplification on both 
the venom gland cDNA of C. squamulosus and F. ornata, as well as that of key outgroup species such 
as Micropechis ikaheka. 

This study reinforces how little we know about the venoms of Australian elapids. Our major 
findings are of interest to the fields of whole organism evolutionary biology and protein evolution, and 
also highlight the hitherto untapped bioresource for drug design and development that the venoms of 
these neglected snakes represent. These findings also have direct implications for treatment of 
envenomed patients. Collectively, these findings strongly support the notion that the venoms of elapid 
snakes underwent a punctuated molecular evolution paralleling the explosive radiation of the snakes 
themselves subsequent to the colonisation of the largely snake free Australian continent 25 million 
years ago. We still have much to learn about the proteomic composition of the venoms of these snakes, 
the role of individual venom components in prey capture and the in vitro activities of these venom 
components. Furthermore, the toxin types recovered in this study should not be considered as the 
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full-suite as this sampling may have only recovered the dominant forms. This study therefore provides 
baseline data for continuing investigations. More extensive sampling is likely to recover novel 
isoforms of toxin types identified to date as well as entirely new toxin classes. In addition, 
investigation of the relationship of the venom gland transcriptomes to the venom proteomes of these 
snakes may reveal translational variances similar to those that have been previously documented for 
other snakes 
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